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Methods for predicting the probability and timing of a species' extinction are typically based on a combination of 
theoretical models and empirical data, and focus on single species population dynamics. Of course, species also interact 
with each other, forming more or less complex networks of interactions. Models to assess extinction risk often lack 
explicit incorporation of these interspecific interactions. We study a birth and death process in which the death rate 
includes an effect from predation. This predation rate is included via a general nonlinear expression for the functional 
response of predation to prey density. We investigate the effects of the foraging parameters (e.g. attack rate and handling 
time) on the mean time to extinction. 

Mean time to extinction varies by orders of magnitude when we alter the foraging parameters, even when we exclude 
the effects of these parameters on the equilibrium population size. In particular we observe an exponential dependence of 
the mean time to extinction on handling time. These findings clearly show that accounting for the nature of interspecific 
interactions is likely to be critically important when estimating extinction risk. 
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1. Introduction 

Population ecologists have long sought to understand 
the relationship between the probability that a popula- 
tion will go extinct, exogenous factors (e.g. environmental 
variability) and endogenous factors (e.g. body size, life 
history, trophic position) (Lande and Stcinar, 2004). As 
well as providing information useful for planning conser- 
vation of individual species, understanding these relation- 
ships promises insights into the determinants of ecosystem 
stability and complexity (McCann, 2000; Allesina et al., 
2012). 

Populations can take different routes to extinction. Ex- 
tinction can be observed as a gradually declining pop- 
ulation size due to habitat deterioration or as a crash 
due to chance events like random catastrophes or random 
fluctuations in population size (demographic stochastic- 
ity) (Lande et al., 1993). Demographic stochasticity is 
caused by random variation in individual mortality and 
reproduction and is independent from time. Environmen- 
tal stochasticity is, on the other hand, random variation in 
time of individual rates (Lande and Steinar, 2004). While 
environmental stochasticity can be important for popula- 
tions of any size, demographic stochasticity becomes par- 
ticularly important at low population sizes. 

Classical population theory shows that demographic 
stochasticity causes the mean time to extinction to in- 
crease exponentially with equilibrium population size, while 
variation in environmental conditions can lead to a power 



law scaling (Lande et al., 1993; Foley, 1994). Ecologists 
have shown which biological features affect species' risk 
of extinction. Factors such as slow life history and small 
geographical range size are all independently associated 
with a high extinction risk (Purvis et al., 2000). Other ef- 
fects such as genetic and environmental change have been 
included into this species-centric approach (Purvis, 2008; 
Lee et al., 2011). One outcome of this large body of re- 
search is population viability analysis (PVA) (Brook et al. , 
2000; Beissinger, 2002). A related outcome is the classi- 
fication of species into extinction risk categories (Mace et 
al., 2008). However single species models used for assessing 
population's viability often lack explicit incorporation of 
direct trophic interactions(Sabo, 2007; Sabo et al., 2008). 

Interspecific interactions have been widely studied the- 
oretically and experimentally in the fields of population 
and community ecology. The pioneering work of Holling 
(Holling, 1959) led to a non-linear density dependent in- 
teraction term called the Functional Response of Preda- 
tion and the effect of functional response on the stabil- 
ity of communities has been widely studied (Drossel et 
al., 2004). There have been several modifications to the 
original Holling's functional response, describing in detail 
different foraging mechanisms (Real, 1977; Abrams et al., 
2000; Jeschke et al., 2002). It has been shown how the in- 
troduction of a scaling exponent between predator attack 
rate and prey density can drastically alter the dynamics of 
realistic systems. This scaling exponent stabilizes chaotic 
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dynamics and reduces or eliminates extinctions in deter- 
ministic food web models (Williams et al., 2004). 

Interspecific interactions could also influence extinction 
time. For example, predation pressure could reduce pop- 
ulation size and thereby increase the chance of extinction 
via variation in population size resulting from demographic 
stochasticity. Predators with different foraging behaviors 
may reduce or increase the risk of extinction of their prey 
according to their feeding strategies. Most PVAs assume 
that any important effects of interspecific interactions are 
summarized in population level parameters (Sabo, 2007), 
for example, predation is included as a constant (den- 
sity independent) source of mortality rather than a cou- 
pled, density dependent population process (Sabo et al., 
2008). Single species models used in PVA often fail to sep- 
arate stochastic variation from population cycles induced 
by species interaction, and this can bias forecasts of via- 
bility (Sabo, 2007; Sabo et al., 2008). 

In this paper we present a single species model in which 
the demographic effects of predation on small prey popula- 
tions are explicitly investigated. We aim to clarify the ef- 
fects of interspecific interactions on the extinction process. 
We analyze a single species birth-death process in which 
the death rate includes density dependent mortality in the 
form of predation by a fixed abundance of predators. We 
then investigate the importance of the predation functional 
response parameters on the mean time to extinction. This 
investigation is novel in itself, and also represents an initial 
step towards a more complete appreciation of the effects 
of interspecific interactions on time to extinction. We use 
the single species framework in order to gain insight into 
the behavior of more complex models. 

2. Stochastic and deterministic models 

We describe the dynamics of single species population 
with number of individuals n(t) at time t by defining func- 
tions for birth rate b(n) and death rate d(n) as 



interaction (Jeschke et al., 2002). We used the general ex- 
pression of the type III functional response 



b(n) 
d{n) 



n 
~ k 
f(n)p, 



(1) 



The logistic form of the birth rate function represents den- 
sity dependent effects such as intraspecific competition for 
resources. The parameter A is the per capita birth rate of 
the prey species and k represents the maximum number 
of individuals the environment can support (Nisbet and 
Gurney, 1982; Nasell, 2001). The death rate includes a 
constant term ^i, the predator-free per capita death rate, 
and the term f(n)p, the per capita functional response of 
predation rate to prey abundance. The function p(t) is 
the number of predator individuals at time t. The per 
capita growth rate of the population is then the difference 
between per capita birth and death rates i.e., r = A — fJ,. 

The functional response f(n) specifies the effect of prey 
density on the predation rate. Historically several func- 
tional responses have been proposed to describe this trophic 



f(n) 



1 + hani+ 1 ' 



(2) 



where a is the attack rate (a measure of the encounter rate 
and capture success of the predator foraging on the prey) , 
h is the handling time (a measure of the time needed to 
attack, eat and digest the prey) and q is a scaling exponent 
between encounter rate and prey density. Different values 
of q are generally assumed to represent different types of 
predation. The type III functional response has been as- 
sociated with learning effects of the predator in catching 
and handling its prey (Real, 1977), with generalist preda- 
tors switching among alternate prey (Smout et al., 2010) 
or with spatial effects enabling prey to hide from predators 
(Vucic-Pestic et al., 2010). Setting q = into expression 
(2) we obtain the type II functional response and setting 
q = and h = we obtain the type I functional response 
(Holling, 1959; Real, 1977). 

The range of variation of the foraging parameters is 
constrained by biological arguments: 

• Handling time h takes values between O.OOldays and 
O.ldays. We assume that h cannot be larger than the 
average lifetime of the prey species (assumed to be 
1/A). 

• Attack rate a takes values between Q.Q0lday~ 1 ind~ 1 
and 10day~ 1 ind~ 1 . This choice is justified by empir- 
ical observations of foraging behavior of insects and 
other organisms (Vucic-Pestic et al., 2010; Hammill 
et al., 2010). 

• The exponent q takes values between and 2 again 
due to empirical observations (Vucic-Pestic et al., 
2010) and theoretical work (Williams et al., 2004). 

Combining expressions (1) and (2) we define the pop- 
ulation birth and death rates B(n) and D(n) as 



B(n) = nb(n) = nX (l - |) 



D(n) = nd(n) = nfi + 



apn 



9+1 



1 + ahni+ 1 ' 



(3) 



The state of the system is completely characterized by the 
probability p(n, t) of having n individuals at time t, where 
n takes integer values in the range {0, • • • , k}. The mas- 
ter equation describing the time evolution of this resulting 
probability distribution is 

d ^ t ± =D(n + l)p(n + 1, t) + B(n - l)p(n - 1, t) 
-(B(n) + D(n))p(n,t). 

In order to express equation (4) in a compact way we in- 
troduce the conditions p(k + 1, t) = p(— 1, t) = 0. Once we 
know the initial condition po (n) = p(n, 0) , then the state 
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probabilities are given by the solution of equation (4) . The 
process represented in this master equation has a degen- 
erate stationary distribution p = (1, 0, • • • ,0), which the 
distribution p(n, t) approaches as time t approaches infin- 
ity. In other words, extinction is ultimately a certainty. 
However, in section 2.2 we will show how to compute the 
conditioned probability of having a population of size n 
at time t conditioned on the fact that it has not yet gone 
extinct. In order to understand the qualitative behavior 
of the model it also helps to understand the dynamics of 
the associated rate equation. 

2.1. Deterministic behavior 

Every stochastic birth and death process has an as- 
sociated deterministic rate equation, describing the time 
evolution of the mean population size (Gardinicr, 2009), 
given by 



dn 
~dt 



B(n) - D(n). 



(5) 



In our case, substituting the rates (3) into equation (5) we 
obtain 

dn ( n\ apn q+1 

_ = An ^l__j_ Mn _ i + aftnffi . ( 6 ) 

At this point we assume constant presence of predators 
and set the number of predators p as a constant. This 
choice enable us to investigate the effect of a generalist 
predator whose density is not affected by the prey den- 
sity. In order to understand how the prey population goes 
extinct, in appendix A we derive an adimensional formu- 
lation of equation (6) and compute all its fixed points and 
their stability. 

In the deterministic rate equation (6) the extinction 
state n = is always a fixed point, while the other fixed 
points are given by the intercepts between the per capita 
growth term (r — Xn/k) and the per capita functional re- 
sponse pf(n) (see figure 1). Two different extinction sce- 
narios are revealed by the analysis of the deterministic 
equation, associated with changes in the stability of the 
extinction state (Assaf et al., 2010). When q = (type 
II functional response) the model has at most two non 
zero fixed points and the extinction state can be stable or 
unstable. When q > (type III functional response) the 
extinction state will always be unstable. 

For the type II functional response the extinction state 
is stable if pa > r. In this case there are two positive fixed 
points, one stable and one unstable (see figure 1). At pa — 
r there is a transcritical bifurcation, at which the unstable 
fixed point goes to (see figure 1) and the extinction state 
becomes unstable. When the extinction state is unstable 
(pa < r) there is only one other fixed point and this is 
stable (see figure 1). For the type III functional response 
(q > 0) the extinction state will always be unstable like 
in the case of the type II functional response after the 
transcritical bifurcation. But, differently from the type II 
functional response, for q > there can be up to three 



non zero fixed points (see appendix A). In this case there 
are combinations of foraging parameters which lead to two 
stable fixed points (see figure 1). 

2.2. Stochastic behavior 

Given the birth and death process (equation (3)), ex- 
tinction is caused by demographic stochasticity and typ- 
ically occurs in two different ways depending on whether 
the extinction state is a stable or unstable fixed point of 
the deterministic rate equation (6). For the type II func- 
tional response, before the bifurcation, extinction is caused 
by a large fluctuation which brings the system from the 
stable to below the unstable fixed point. From there the 
fast deterministic evolution takes the population quickly to 
the stable extinction state. For the type II functional re- 
sponse after the bifurcation, and for the type III functional 
response, extinction is caused by a large fluctuation in den- 
sity which brings the system from the stable fixed point 
to the unstable, absorbing extinction state. For the type 
III functional response, when there are two stable fixed 
points, extinction needs two large fluctuations to occur, 
one fluctuation which brings the system from the larger to 
the lower stable fixed point and another fluctuation which 
brings the system from the lower stable fixed point to the 
unstable extinction state. This result confirms the stabi- 
lizing effect of the type III functional response (Williams 
et al., 2004). 

When the deterministic rate equation has at least one 
stable fixed point, the system approaches a quasi station- 
ary state with a time independent distribution 7r(n); this 
is called the Quasistationary Distribution (Bartlctt, 1960; 
Nisbet and Gurney, 1982). The quasistationary distribu- 
tion 7r(n) is obtained from the probability p c (n, t) of find- 
ing n individuals at time t, conditioned on the fact extinc- 
tion has not occurred yet: 



Pc(n,t) 



p(n, t) 
1-P(0,t)' 



(7) 



We derive a master equation for the conditioned probabil- 
ity p c (n,t) and look for its stationary solution n(n) (see 
appendix B). When the initial condition of equation (4) is 
set to the quasistationary distribution, then the probabil- 
ity of finding n individuals at time t becomes 



p(n,t) ~ x(ri)exp(-t/MTE). 



(8) 



The time to extinction is then an exponentially distributed 
random variable with mean equal to MTE and, if the 
initial condition of (4) is set equal to 7r(n), then the mean 
time to extinction is 



MTE 



£>(1)tt(1) 



(9) 



There are more complicated expressions for the mean time 
to extinction when the initial condition is not the quasi 
stationary distribution ir(n) (see appendix B). 
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Figure 1: Prey population mortality due to predation (left panels) and per capita mortality due to predation (right panels) following a type 
II (upper panels) and a type III (lower panels) functional response. On the left plots the continuous lines A, B and C represent different 
population functional responses. On the right plots solid lines are the per capita functional response while the dashed line is the per capita 
growth curve (A = 1.5 day -1 ; fi = 0.5 day^ 1 and k = 150ind). The number of predators is fixed at p = lind. The foraging parameters 
are A: h = O.QOldays; a = 0.5 day~ 1 ind~ 1 . B: h = 0.02days; a = 1 day~ 1 ind~ 1 at the transcritical bifurcation. C: h = 0.03days; 
a = 2.05 day ~ 1 ind~ 1 . Note that with this particular choice of foraging parameters the value of the fixed point is fixed at nO = bOind. 
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We compute the MTE of the birth and death pro- 
cess (3) for different functional responses and for different 
values of the foraging parameters within the functional re- 
sponses. All other parameters remain fixed. 

2.3. Numerical calculations 

Despite its apparent simplicity it is not possible to ob- 
tain closed expressions for the quasistationary distribution 
n(n) of the birth and death process (3). Instead we ob- 
tain the quasistationary distribution in a realistic range 
of foraging parameters by an iterative numerical scheme 
described in appendix B. In order to perform numerical 
calculations we fix the growth parameters of equation (6) 
in the following way: A = 1.5 day" 1 ; fi — O.dday" 1 ; 
k = 150 ind. With that choice the intrinsic growth rate 
of the prey population is fixed to lday" 1 . If we look at 
equation (6) without predators i.e., putting p — 0, we ob- 
tain a modified version of the classic logistic equation: 



dn 
7/7 



Xn 2 

IT' 



(10) 



The modified logistic equation (10) has a non zero fixed 
point n* — kr/X — 100 ind. The system behaves according 
to equation (10) when predation is absent or very weak. 

We fix the number of predators to p — 1 ind and com- 
pute the logarithm of the Mean Time to Extinction (MTE) 
of the birth and death process (3) obtained using the quasi 
stationary distribution as initial condition (see appendix 
B) for different combinations of foraging parameters. In 
the case of the type II functional response (q = 0), we 
also compute the MTE using different choices of initial 
conditions. In this case the time needed to relax to the 
quasi stationary distribution is negligible with respect to 
the MTE. 

In order to quantify the difference between the effects 
of the foraging parameters on the equilibrium population 
density and the direct effect of the foraging parameters 
on the mean time to extinction we also perform numerical 
simulations of the model changing the foraging parameters 
and keeping fixed the stable fixed point. As we increase 
the handling time, we also increase the attack rate in or- 
der to keep the effective predation pressure constant at 
the equilibrium population. Using the equilibrium rela- 
tion B(n ) = D(n ), we find a relation between attack 
rate and handling time which includes one fixed point no 
as another parameter 



a{h, no) 



rk — Xno 



n^[pk — hno(rk — Xno)] 



(11) 



In order to avoid fixing an unstable equilibrium with re- 
lation (11), we limit our investigation to those values of 
handling time which give raise to a negative Jacobian (see 
appendix A) i.e., 



h < hi = 



Xpk[(l — q)Xno + qrk] 
n (l + q)(rk - Xn ) 2 ' 




log(MTE) 



0.02 0.04 0.06 
Handling Time (days) 



Figure 2: Variation of the stable fixed point no predicted by equation 
(6) (upper panel) and the log of the mean time to extinction (MTE) 
of the birth and death process (1) using the quasistationary distribu- 
tion as initial condition (lower panel) as a function of handling time 
and attack rate for type II functional response. The growth param- 
eters and the number of predators are as specified in the legend of 
figure 1. 



As the handling time approaches the value 

pk 



h 



no(rk — Xno) 



(13) 



(12) 



the required attack rate approaches infinity, so we limit 
ourselves to h < min(ho, hi). In the case of q > we also 
keep h small enough so that we do not enter the bistable 
region. Matlab code to reproduce our calculations can be 
found at http : //purl . org/net/extinction_code. 

3. Results 

The foraging parameters of the functional response af- 
fect the MTE in direct and indirect ways. The foraging 
parameters affect the MTE indirectly by influencing the 
equilibrium population size of the prey. They influence 
it directly by changing the MTE even for a fixed equi- 
librium prey population size. The indirect effects of the 
foraging parameters on MTE can be seen in figure 2. As 
expected, the MTE is relatively low when the functional 
response parameters lead to extinction being predicted by 
the deterministic model (region labeled no = in figure 
2). When the foraging parameters are such that there is a 
stable positive equilibrium, the MTE is clearly positively, 
but non linearly, related to the equilibrium population size. 
Changes in both attack rate and handling time strongly 
influence MTE through changes in equilibrium population 
size. 

Figures 3 and 4 show the direct effect of the foraging 
parameters on the MTE when the equilibrium density is 



kept fixed. In that case we see that, when the extinc- 
tion state is unstable, the MTE decreases exponentially 
with handling time (figure 3), or, in the case of type III 
functional response, more than exponentially (figure 4). 
When the extinction state is stable, the MTE decreases 
less than exponentially (figure 3). We fitted exponential 
curves to the MTE for four different values of n , with type 
II functional response, before the transcritical bifurcation 
i.e. when the extinction state is unstable. We observe an 
increase in both the slope and the intercept with increasing 
carrying capacity 1 . 

The MTE of the model without predation is extremely 
large (10 25 days), meaning that extinction will never oc- 
cur in a biological time when there are no predators. We 
observe a decrease of 20 orders of magnitude in the MTE 
caused by predation interaction. This decrease is partly 
as a result of a decrease in the equilibrium population size 
alone, making extinction due to demographic stochasticity 
more likely. However demographic stochasticity is impor- 
tant also when the equilibrium density is kept fixed. For 
example, for a realistic range of foraging parameters, for 
a type II functional response, at a given fixed point (e.g. 
n = 50) the MTE varies by about 5 orders of magnitude 
with handling time. The magnitude of this variation in- 
creases up to 10 orders of magnitude for increasing values 
of the stable fixed point (figure 3). 

4. Discussion 

We have shown how different choices of the foraging 
parameters vary the mean time to extinction by up to 10 
orders of magnitude when equilibrium population is kept 
constant. The effect of the foraging parameters can be 
intuitively explained by looking at the strength of popula- 
tion regulation at the equilibrium density i.e., the slope of 
the functional response at the fixed point (Figure 1). Dif- 
ferent handling times and attack rates that keep the stable 
equilibrium fixed result in decreasing maximum saturation 
rate (increasing handling time) and increasing half satu- 
ration density (increasing attack rate). This regulation 
doesn't alter the actual carrying capacity of the prey pop- 
ulation, but can change drastically the demographic effects 
of predation. These demographic effects are independent 
of environmental stochasticity and become relevant for ex- 
tinction risk only for small populations sizes (Lande and 
Steinar, 2004). Therefore when the predation interaction 
is present it is not possible to reliably estimate extinction 
risk without having a measure of the foraging parameters. 

There is a wide literature describing experimental mea- 
sures of foraging parameters (attack rate, handling time 
and scaling exponent) in laboratory communities. These 
studies include predator prey interactions among terres- 
trial and aquatic organisms such as protists (Hammill et 



1 The values of the slope are fitted using a least square method 
and are, for different values of the fixed point: no = 40, slope -85.8; 
no = 50, slope -133.7; no = 60, slope -168.3; no = 70 slope -173. 



al, 2010) and arthropods (Spitze, 1985; Smout et al., 2010). 
However similar measurements of the nature of predator 
prey interactions are absent in most of the studies related 
to the extinction risk of individual species. The previous 
focus on individual species characteristics is presumably 
due to lack of multispecies time series data (Sabo et al., 
2008). 

We have shown with a simple model how trophic in- 
teraction can be relevant in assessing extinction risk of a 
target species. Our result could be investigated analyti- 
cally using refined approximation techniques (Ovaskaincn 
et al., 2010). In (Nasell, 2001) and (Assaf et al., 2010) and 
an approximate expression of the quasistationary distri- 
bution and the mean time to extinction is derived for the 
stochastic logistic model and the SIS model of epidemics, 
both slight simplification of our model. The bistability 
emerging for high values of handling time, for type III 
functional response, could also be investigated using ap- 
proximation techniques. 

Another natural route for further investigations would 
be to constrain our analysis of parameter space to combi- 
nations of foraging parameters that occur in reality. Intro- 
ducing allometric relationships between foraging parame- 
ters and relating them to the growth parameters would be 
one way to constrain such analysis. We have shown this 
dependence as a function of the stable equilibrium density 
in equation (11). This relation can be generalized using 
allometric scaling relations between attack rate and han- 
dling time (Brose et al., 2006). Such an allometric scaling 
would require a more general formulation of the model, 
including predator biomass and prey biomass as other pa- 
rameters. 

This work has application both on studies on extinc- 
tion risk and on studies on foraging theory and we aim 
to extend our analysis to multispecies communities in fu- 
ture studies. Most of the existing theoretical studies about 
complex communities do not incorporate the effects of de- 
mographic stochasticity and use deterministic measures of 
persistence to assess the extinction risk (Hofbauer et al., 
2008; Dunne ct al., 2009; Sahasrabudhe et al., 2011). Our 
findings fit into recent application of stochastic methods 
to population biology (McKane et al., 2004; Ovaskainen et 
al., 2010; Black et al., 2012) and can be used to increase 
the predictive understanding of extinction processes. 

For further information and discussion of this paper see 
http : //purl . org/net/ extinction. 

A. Qualitative Behavior 

We derive an adimensional formulation of (6) . We scale 
the number of individuals with the actual carrying capac- 
ity m = n/k, and the characteristic time with the intrinsic 
birth rate t = At. With this substitution we have: 

drn(r) m q+1 

— = mil - m) : —r. (14) 

dr v ; a + bmi +1 y ' 
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Figure 3: Logarithm of the mean time to extinction (MTE) (when the quasistationary distribution is set as initial condition of equation (4)) 
as a function of handling time for different values of the stable fixed point no, for a type II functional response. The dashed line shows the 
values of handling time at which there is a change in the stability of the extinction state (left is unstable, right is stable). The curves are 
drawn for all values of handling time which keep the stable equilibrium density no fixed and the corresponding attack rate (11) finite i.e., 
for h < min(hi,ha). The values of /ii and /io are obtained from expressions (12) and (13). We used expression (9) to obtain the MTE. The 
growth parameters and the number of predators are as specified in the legend of figure 1. 




Figure 4: Logarithm of the mean time to extinction (MTE) as a function of handling time for different values of the stable fixed point no, for 
type III functional response (q = 1). We used expression (9) to obtain the MTE. The growth parameters and the number of predators are as 
specified in the legend of figure 1. 
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Extinction is always a fixed point of the system. The 
other fixed points are where the per capita growth term 
equals the per capita predation rate and are given by the 
intercepts between the per capita growth term (J — n) and 
the per capita functional response [m q /(a + bm q+1 )\ (see 
figure 1). In mathematical terms, putting dm/dr = 
into equation 14, the fixed points of the system are the 
extinction state n = and the solutions of the algebraic 
equation 



bm q+2 - blm q+1 



la = 0. 



(15) 



Where the adimensional parameters have following mean- 
ings: 



• 7=1 — l/i?o where Rq 
ductive ratio 



A//i is called basic repro- 



• a is the ratio between intrinsic birth rate and per 
capita attack rate at the maximum population den- 
sity of the prey species 

A 



a(q) 



ak q p 



(16) 



Note that a is the only adimensional parameter which 
depends on the scaling exponent q. This parameter 
can be called inverse adimensional attack rate. 

• The last parameter is the one which takes into ac- 
count handling time h 

k 



Xh- 



P 



(17) 



This can be seen as the ratio between handling time 
and average generation time scaled with the ratio 
between maximum prey density k and the (fixed) 
number of predators p. This parameter can be called 
scaled adimensional handling time. 

From the right lower panel of figure 1 we see that equa- 
tion (15) has at most 3 real and positive solutions. Note 
that there will be positive solutions only when < Rq < 1 
i.e. when A < \i. 

As a special case we can look at the solutions of (15) 
when <7 = 0, i.e., for type II functional response. Then 
equation (15) simplies to 

bm 2 + (a- bl)m+ (1 - al) = 0. (18) 

This equation has real solutions when (a + bl) 2 > 46. If 
this is true then the two real solutions of (18) will be 



mi 9 = 



(bl -a)± y/(a + bl) 2 
2(1 - al) 



46 



(19) 



Now we can look at the stability of the fixed points of 
equation (14), performing the classical stability analysis. 
We compute the Jacobian of equation (14) 

d ( dm\ (q + l)am q 

dm 



J(m) 



dt J {a + bm q + 1 ) 21 

and putting m — into (20) we have 



(20) 



• if q — then J(0) = r — ap where r = A — [i so 

1. if r < ap the extinction state is stable 

2. if r > ap the extinction state is unstable 

• if q > the extinction state is always unstable 

The stability of the extinction state gives raise to the the 
two extinction scenarios described in the Methods section. 

B. Mean time to Extinction 

In this appendix we present the mathematical and nu- 
merical tools needed to obtain the quasi-stationary distri- 
bution and the mean time to extinction of a general birth 
and death process when ultimate extinction is certain. We 
define the process as the time evolution of a random vari- 
able {X(t), t > 0} in a finite space state {0, 1, • • • , k} where 
the origin is an absorbing barrier (Nisbet and Gurney, 
1982). We recall the master equation (4) and write it in a 
more compact way using a vectorial notation: 



dp(t) 
dt 



p(t)A, 



(21) 



where p(t) = (p(0, t),p(l, t), ■ ■ ■ ,p(k,t)) is the row vector 
containing the state probabilities and the matrix A con- 
tains the transition rates as follows: 



A = 



/ -G(0) 
D(l) 


V o 



B(0) 
-G(l) 
D(2) 







5(1) 
-G(2) 





(22) 



-G(k) J 



with G(n) = B(n) + D{n). A is a tridiagonal matrix with 
all row sums equals to 0. Note also that, using the rates 
(3) the first row is a row of zeros. The solution of the 
master equation (4) will give the probability of having n 
individuals at time t; in other words p(n,i) = P{X{t) = 
n}. Once the birth and death process is defined we de- 
rive two auxiliary processes (Nasell, 2001) {X^°'(t}} and 
{X^>(t)} associated to {X(t)}, both with reduced state 
space {1, 2, • • • , fc}: 

• The first process {I<°)(t)} is equal to {X{t)} but 
with reduced state space. So the rates of transition 
are not changed except D(l): 



B^\n) = B(n), 
D(°)(n) = D(n), 
= 0. 



(23) 



But the assumption on D^(l) ensures that there 
will never be extinction, in other words that there is 
no absorbing state. 

• The second process {X^\t)} is defined, in the re- 
duced state space, with the following transition rates: 



B^(n) = B(n), 
D^(n) = D(n-1), 
£>«(0)=0, 



(24) 



and corresponds to a population in which there is 
always at least one "immortal" individual. 

It is important to underline that the two auxiliary pro- 
cesses we have defined have no ecological meaning. {X^ (t)} 
and {X^(t)} will be used only as mathematical tools in 
order to obtain information on the time to extinction of 
the birth-death process (1). 

The stationary distributions of the two auxiliary pro- 
cesses are: 



pW(n) 



T(n) 



(25) 



where T(n) and R(n) are defined as follows: 



T(n) 
R(n) 



B(1)B{2) ■ ■ ■ B{n - 1) 
D{2)D{Z) ■ ■ ■ D(n) ' 
B{l)B(2) ■ ■ ■ B(n - 1) 
D(l)D{2) ■ ■ ■ D{n - 1)' 



(26) 



Note that T(n) — j <-w D(n) • 

Now we can partition the state space of the original 
process into two subsets, {0} and Q = {1, 2, • • ■ , fc}. Q is 
the set o transients for {X(t)} and the state space for the 
two auxiliary processes while {0} is the absorbing state 
for {X(t)}. Correspondingly we can partition the state 
vector p(f) and the transition matrix A, and obtain from 
equation (21) 



dp(0,t) _ dp Q (t) 



dt 



dt 



[p(0,t); PQ (t)] 







(27) 



Here Pq(£) is the vector of probabilities in the transient 
states and a = (D(l), 0, • ■ ■ ,0). With this separation we 
can split the master equation (21) into: 



dp(0,t)/dt = p Q (t)a = D(l)p(l, t), 
dp Q (t)/dt = p Q (t)A Q . 



(28) 



Before absorption the process takes values in the set of the 
transients. 

As in equation (7) we define the conditional probability 
p c (n,t) = P{X(t) = n\X(t) > 0} of having n individuals 
at time t knowing absorption has not occurred yet. And, 
using equations (28), we can express it in vectorial form: 



Pc(t) 



i-p(o,i) i- P (p,ty 



p(t) 



(29) 



Differentiating equation (29) and using the master equa- 
tion (4) and equations (28) we obtain an equation for p c (t): 



dp c (t) _ dp Q ( 1 



p Q {t) dp(fl,t) 



dt dt \l-p(0,t)J (l-p(0,i)) 2 dt 

= p c (t)A Q + D{l) Pc (l,t)p c (t). 



(30) 



Putting to zero the right-hand side of expression (30) we 
obtain an equation for the quasi-stationary distribution 
7r = (tt(1), 7r(2), • • • ,7r(fc)), defined as the distribution of 
the transient states conditioned on the fact that there has 
been no extinction yet: 



ttAq = -D(l)7r(l)7r. 



(31) 



In other words the quasistationary distribution it is the 
left eigenvector of Aq with eigenvalue — D(1)tt(1). 

It can be shown (Nasell, 2001) using a recursive solu- 
tion of the master equation (30), that the n(n) satisfy the 
recursive formula: 



n{n) = T(n) ^ 



(i-E^U(j)) 
R(i) 



tt(I). 



(32) 



Once 7r(l) is known then 7r(2), 7r(3), . . . , ir(n) can be deter- 
mined iteratively. But 7r(l) can only be obtained knowing 
all the other elements by the relation E ra 7r ( n ) = 1- ^ or 
this reason the analytic determination of 7r is limited to 
birth death processes with linear transition rates and this 
is not our case. However there is an iterative method that 
can be used to derive numerical approximations for the 
quasistationary distribution of our process: 

• Start with an initial guess for 7r(l). 

• Obtain all the 7t(n) using the (32) and compute S = 

• Start the iteration again with 71^(1) = ir(l)/S and 
obtain all the 7r / (n). 

• Repeat the process until \\ir K+1 (n) — 7r^(n)|| < 5. 
The value 5 gives the precision of the algorithm. 

The time to extinction 2 r is a random variable that 
depends on the initial distribution of the process (Nasell, 
2001). We call tq the time to extinction of the birth death 
process when the quasi stationary distribution 7T is set as 
initial distribution. And we call r„ the time to extinction 
when the initial condition is X(0) = n i.e., when p(n, 0) = 
1. 

If absorption has occurred at time t then the events 
{t < t} and {X(t) — 0} are identical: 



P{T<t} = P{X(t)=0}= P (0,t). 



(33) 



Once the quasistationary distribution is known then the 
MTE 3 is given by expression (9). 

The explicit expression of time to extinction with an 
arbitrary initial condition is more difficult to obtain. It is 
a standard result for birth death processes theory (Nisbet 



2 Be careful this r is different from the adimensional time used in 
appendix A 



'note that with our notation MTE = E(tq). 
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and Gurney, 1982) that at least the expectation can be 
determined explicitly when X(0) = n: 

-j n k 
^n) = W Y. W) Y.TW- (34) 

Note that putting n — 1 in the above formula the expected 
time to extinction from the state 1 can be written as fol- 
lows: 

k 

E ^ = W)^ m = w^ry (35) 

The expected time to extinction from a state n can there- 
fore be written in an alternative form in function of the 
stationary distribution of the first auxiliary process: 

n k 

Then the expected time to extinction for an arbitrary ini- 
tial distribution {p(n, 0)} can be derived from (34): 

^) = 1 ^ ) tm± i ^± P (n,0), (37) 

j — 1 i— 1 n— i 

with the assumption that the initial distribution is sup- 
ported on the set of the transient states. 
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